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Abstract 

In a recent paper jTj we have presented an automated subtraction method 
for divergent multi-loop /leg integrals in dimensional regularisation which allows 
for their numerical evaluation, and applied it to diagrams with massless internal 
lines. Here we show how to extend this algorithm to Feynman diagrams with 
massive propagators and arbitrary propagator powers. As applications, we 
present numerical results for the master 2-loop 4-point topologies with massive 
internal lines occurring in Bhabha scattering at two loops, and for the master 
integrals of planar and non-planar massless double box graphs with two off- 
shell legs. We also evaluate numerically some two-point functions up to 5 
loops relevant for beta-function calculations, and a 3-loop 4-point function, the 
massless on-shell planar triple box. Whereas the 4-point functions are evaluated 
in non-physical kinematic regions, the results for the propagator functions are 
valid for arbitrary kinematics. 



1 Introduction 



The precision of present and upcoming high energy experiments is unprecedented 
and clearly requires an analogous level of accuracy on the theoretical side. One of the 
steps towards achieving this goal is to push the level of higher order corrections in 
perturbation theory further. However, the calculation of the corresponding Feynman 
diagrams becomes more and more involved as the number of loops and legs of the 
graphs increases. Typically the occurring tensor integrals are expressed in terms of 
scalar integrals which are then reduced to a set of master topologies. These master 
topologies are the basic ingredients and thus the bottleneck for any attempt to perform 
a higher order calculation. 

A milestone in this respect has been the analytic calculation of the master integrals 
for the planar |2j and non-planar j ! 3 j massless two-loop box diagrams with on-shell 
external legs and the ones with one external leg off-shell [U EJ El H] • These results 
triggered the analytical computation of two-loop matrix elements for 2^2 scattering 
processes with on-shell legs [S] and 1 — > 3 processes with one off-shell leg j^j, which are 
necessary for the next-to-next-to-leading order calculation of prominent processes like 
two-jet production in hadronic collisions or three-jet production in e + e~ collisions. 
On the other hand, the two-loop four-point functions with two external legs off-shell 
still await their analytical calculation. These integrals are for example needed to 
calculate the production of two massive vector bosons in hadronic collisions at two 
loops. 

The reason for the lack of analytical results is the increasing number of scales as 
the number of off-shell legs increases, leading to more complicated analytic structures. 
The same argument is true if the diagrams contain massive internal lines. Such dia- 
grams are certainly important in view of a forthcoming e + e~ collider at TeV energies, 
where the experimental precision will be such that the Standard Model is tested at the 
two-loop level. However, a major part of the corresponding two-loop master integrals 
with massive propagators may not be accessible analytically. Thus the development 
of numerical methods seems natural to make progress in loop calculations. 

This is not at all straightforward, as different kinds of singularities - infrared 
(IR), ultraviolet (UV) and threshold singularities - are present. Even in the absence 
of IR and UV divergences, one has to face problems due to a complicated analytic 
structure: For physical kinematics the integral generically has singularities within the 
integration region, which typically hinder a successful, numerically stable evaluation. 
At the one-loop level this problem is less severe because the integration space can be 
reduced to lower dimensionality and because the analytical structure is completely 
understood. Solutions for multi-leg one-loop diagrams have been suggested in [TUJUI]. 
However, beyond one loop the problem is not satisfactorily solved in a general and 
constructive way, although promising approaches have been proposed [121 EB1 EH HHj • 
In the presence of infrared divergences the situation is even more difficult. One first 
has to obtain a numerically integrable function by subtracting the infrared poles. As 
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these poles are in general overlapping, this is a highly non-trivial task. In this paper 
we will focus on this problem, i.e. on the isolation of infrared divergences from general 
multi-loop Feynman diagrams. 

In P] we proposed a constructive method which enables to calculate multi-loop 
integrals numerically for unphysical kinematics, and applied it to diagrams with mass- 
less propagators. It served to check the results for the massless on-shell double boxes 
and to make predictions for the massless doubles boxes with one off-shell leg, the 
latter being confirmed by the subsequent analytical calculation later. This method is 
based on sector decomposition, originally used to disentangle overlapping UV diver- 
gences in the proof of the BPHZ theorem ^Bj. The basic concept being very simple, 
it has a wide range of applicability. For example, it has been used to derive high en- 
ergy approximations for one-loop integrals ^7j. Sector decomposition also has been 
employed to extract logarithmic mass singularities from massive multi-scale integrals 
in the high energy limit at two loops ^Hj- In 0, the concept of sector decomposition 
has been elaborated to a highly automated program package, allowing to convert di- 
mensionally regulated Feynman integrals into a Laurent series in e. The coefficients 
of the poles are sums of multi-dimensional parameter integrals. In a kinematic region 
where thresholds are absent, the latter can be evaluated with standard numerical in- 
tegration routines. In simple cases, analytic integration is also possible. The numbers 
of loops and legs are in principle arbitrary, but of course limited by CPU time and 
disk space. Recently the method has been extended to deal also with phase space 
integrals [19J. 

In this article, we show how to generalise the proposed method to arbitrary Feyn- 
man diagrams. This is achieved by extending the algorithm to graphs with arbitrary 
masses and propagator exponents. With these extensions one is able to deal with 
very different types of multi- loop/multi-scale problems, most of which are at the edge 
or beyond the present state of the art of analytic calculations. More precisely, we 
present numerical results for massless two-loop 4-point functions with two external 
legs off-shell, and for the three most difficult master topologies for Bhabha scattering 
at two loops, at two numerical points where s,t < 0. We also give numerical results 
for 3-loop, 4-loop and 5-loop massless 2-point functions relevant for the calculation 
of beta-functions, and finally for the massless on-shell planar 3-loop box diagram. 
Although we do not provide a solution for arbitrary physical kinematics in general, 
we note that there is an important class of diagrams, namely diagrams with only one 
scale, where there is no kinematical restriction. 

The paper is organised as follows. In section 2, we briefly present the method, 
extended to include also non-integer propagator powers, respectively Feynman param- 
eters in the numerator. In section 3 we collect applications of the method, presenting 
numerical results for the types of diagrams listed above. Section 4 contains the dis- 
cussion of our results and the conclusions. 
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2 The formalism 



The formalism of iterated sector decomposition enables to extract in a constructive 
way the divergences of arbitrary scalar L-loop iV-point integrals, where the diver- 
gences are regularised by dimensional regularisation. The method has been described 
in detail in £[]. Here we only sketch it briefly, focusing on the features that have not 
been elaborated previously, i.e. the possibility of arbitrary propagator powers and 
the inclusion of massive propagators. 

A scalar graph G with N propagators and L D-dimensional loop momenta, where 
the propagators can have arbitrary, not necessarily integer powers i/j , has the following 
representation in momentum and Feynman parameter space 1 : 

f L d D ki N 1 
G = J H jS (q] ' m) + W> 

v(n -rn/^7 N , N jjn v -(l+i)d/2 

iij=ii \y 3 ) J j= i i= i s 

Here qj are the propagator momenta, i.e. linear combinations of external and loop 
momenta, N u = YljLi v j- The functions Li and T can be straightforwardly derived 
from the momentum representation. They also can be constructed from the topology 
of the corresponding Feynman graph as follows. 

Cutting L lines of a given connected L-loop graph such that it becomes a connected 
tree graph T defines a chord C(T) as being the set of lines not belonging to this tree. 
The Feynman parameters associated with each chord define a monomial of degree L. 
The set of all such trees (or 1 -trees) is denoted by 7[. The 1-trees T e T\ define IA as 
being the sum over all monomials corresponding to a chord C(T G %). Cutting one 
more line of a 1-tree leads to two disconnected trees, or a 2-tree T. T 2 is the set of 
all such 2-trees. The corresponding chords define monomials of degree L + 1. Each 
2-tree of a graph corresponds to a cut defined by cutting the lines which connected 
the 2 now disconnected trees in the original graph. The momentum flow through the 
lines of such a cut defines a Lorentz invariant Sf = (X^- 6 cut(T ) Vj) 1 - The function JF 
is the sum over all such monomials times minus the corresponding invariant: 



u{x) = e [ n 

TeTi jeC(T) 

*b(a) = E [ II 

fer 2 j€C(f) 



-Sf) , 



N 

F(x) = JF (x) +U{x)Y J x ] m) . (2) 

3=1 

1 This representation, as well as the following discussion of its topological properties, is well 
known, see for example [221 051 121 1251 IT5] . 
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IA is a positive semi- definite function. Its vanishing is related to the UV sub divergences 
of the graph. Overall UV divergences, if present, will always be contained in the 
prefactor T(N U — LD/2). In the Euclidean region, T is also a positive semi-definite 
function of the Feynman parameters Xj. Its vanishing does not necessarily lead to 
an IR singularity. Only if some of the invariants are zero, for example if some of the 
external momenta are light-like, the vanishing of T may induce an IR divergence. 
Thus it depends on the kinematics and not only on the topology (like in the UV case) 
whether a zero of T leads to a divergence or not. This fact makes it much harder 
to formulate general theorems on the IR singularity structure of Feynman graphs. 
Examples of U and T functions will be given below. 

In multi-loop integrals, the singular regions - leading to at most l/e 2L infrared 
poles upon integration in parameter space - are generally overlapping, such that the 
set-up of a simple subtraction scheme becomes impossible. In pQ we have developed 
an algorithm and computer program to disentangle such overlapping regions. It can 
briefly be sketched as follows: 

Step 1: A "primary" sector decomposition of all Feynman parameters eliminates the 
5-function and maps the integral to a sum of iV (iV-l)-dimensional parameter 
integrals over the unit cube: The decomposition 

N „ N 



Jal N x = J2j al N x J] 0{xi > x j > 0) 



1=1 3=1 

and appropriate variable substitutions lead to 

N IjV-l j,N-(L+l)D/2 



G = ^Gi, G,= / I] dx i 
1=1 o i=1 



^- LD/ \x) 



where lAi,Ti contain only positive semi-definite functions and where the singu- 
larities can only be located at Xj = 0. 

Step 2: Sector decomposition is performed iteratively for sets of parameters {xk} 
which make the transformed T or U functions vanish at x^ = 0. The iteration 
stops if U — 1 + . . . and T = (— Sy) + . . . , where the dots denote positive semi- 
definite remainder terms and (— s^-) is a kinematic invariant which is assumed 
to be positive. At this point all divergences are non-overlapping, the poles will 
be entirely contained in overall factors with a a negative integer. 

Step 3: Subtractions are performed using Taylor expansion around Xj = 0. After 
subtraction, it is safe to expand in e up to the desired order, leading to the form 

G = £ ^M + (e M+1 ). (3) 

j=-M €l 
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Note that the method in principle allows to calculate coefficients of the Laurent 
series in e to arbitrary order M. 

Step 4: The coefficients C m (x) in Q are sums of subtracted and expanded subsector 
integrals in terms of bounded functions of Feynman parameters. These can be 
evaluated numerically, or even analytically in simple cases. 

In we were concentrating on massless integrals. The iteration always stopped, 
as each sector decomposition reduces the degree of the monomials which are relevant 
for the IR behaviour of the integral. Including propagator masses, it is clear that 
the IR behaviour of a diagram can only improve, and hence this extension of the 
algorithm does not pose any principle problem. On the other hand, a new technical 
feature appears, which is that the Feynman parameters associated with massive lines 
occur quadratically in J 7 , stemming from the term proportional to U in eq. (0). In 
this case it is not guaranteed that the iteration stops if one proceeds thoughtlessly. 
For example, it is easy to see that a polynomial of the type xy 2 + z 2 can produce 
an endless loop if a sector decomposition in {x, z} is followed by one in {x,y}. A 
possible way out is to decompose first all parameters which occur quadratically, in 
our example {y, z}. Then no dangerous situation is present anymore and one can 
proceed as usual. The corresponding modification of the code is straightforward, and 
examples of such a situation are the two-loop QED graphs calculated below. 

We note that with the generalisation of the algorithm presented in this article, 
integrals with arbitrary numerators, and thus tensor integrals, can be treated. As 
is well known, after integration over the loop momenta, tensor integrals correspond 
to linear combinations of integrals with monomials of Feynman parameters in the 
numerator (times tensors carrying the Lorentz structure). These numerators are 
taken into account by our program and may cancel IR singularities which would 
be present in the scalar diagram. Another possibility is to express tensor integrals 
by Feynman parameter integrals which do not have polynomial numerators, but are 
partly in shifted dimensions, as advocated in |25l 12*5] . As our method treats the 
space-time dimension as a parameter, this way is viable as well. 

3 Applications 

In this section we discuss a number of applications of our algorithm, ranging from 
2-loop 4-point functions to 5-loop 2-point functions. We note that sector decom- 
position also is useful for the calculation of graphs which are not plagued by IR or 
UV divergences, as it produces integrands which are bounded and positive definite 
in kinematic regions where no thresholds are approached. The technical problem of 
integrable singularities inside multi-dimensional integration regions hinders the direct 
computation of the parameter integrals for general kinematics. However, this is not 
a limitation intrinsic to our procedure, as more sophisticated numerical integration 
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routines than the one we use could overcome this problem. However, for the time 
being, we have to restrict ourselves to the non-physical region, where all Mandelstam 
variables are negative, in the case of the box functions. 

For the numerical evaluations we have used the standard Monte Carlo package 
BASES [SHI- The numerical precision for all quoted results is 1% or better. 

The integration time for all 2-loop boxes to reach the demanded precision is of 
the order of an hour on a PC with a Pentium IV (2 GHz) processor. For the 3-loop 
propagators the required time is less, whereas for the 5-loop propagator it took a few 
days to reach this precision. 

3.1 Parameter representation of general double box graphs 

As the topological functions U and T are the basic input for our algorithm, we first 
give the parameter representation of the general planar and non-planar double box 
topologies. These representations are used as a starting point for the particular 2-loop 
box examples calculated below. 

Planar topology 

The functions U and T for the planar double box graph Gp(s, t, u, s%, s 2 , S3, S4, {m 2 }), 
as shown in Fig. (JT^), are given by 

U — ^123^567 + ^4^123567 

F = (~ S ) (^2^3^4567 + ^5^6^1234 + X 2 X 4 X 6 + £3X4X5) + (— t) £1X4X7 
+ (~Sl) Xi (X 4 X 5 + X 2 X 4567 ) + (-8 2 ) Xi (X 4 X 6 + X 3 X45 67 ) 

+ (-S3) £7 (£3£4 + £e£i234) + {-Si) x 7 (x 2 x A + £5X1234) 



where we use the short-hand notations = p 2 , s = (p\ +p 2 ) 2 ,t = (P2 + P3) 2 and 
x h...i n — £«i + • • • + %i n - In the examples considered in the following, we will set some 
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(4) 



of the Si and/or m 2 to zero. 



Non-planar topology 

The parameter representation of the general non-planar double box 
GWO, t, u, si, s 2 , s 3 , s 4 , {m 2 }), Fig. ^TJb), is given by 



U 



£l23£4567 + £45£67 

(-S) (£ 2 £3^4567 + £2£4£6 + £3^7) 
+ (—t) X1X4X7 + (—u) X1X5X6 



T 
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+ ("Si) £1(25X7 + X 2 X 45 q 7 ) + (~S 2 ) Xi(x a Xq + £3X4567) 
+ (-S3) (x 6 X 7 Xi2345 + + X 2 X 5 X 6 ) 

+ (-S4) (^4^5^12367 + a^a^e + X 2 £ 4 X 7 ) 

7 

+ W^i jm j. (5) 

3.2 Massless double box graphs with two legs off-shell 

For the massless 4-point functions considered in this subsection, two legs are off-shell 
and the internal lines are massless. The topologically different positions of the off- 
shell legs give rise to three master integrals for both the planar- and the non-planar 
double box with two legs off-shell. We give results at two Euclidean numerical points 
for each master integral, as shown in Tables ^ and |21 An overall factor T 2 (l + e) has 
been extracted: 

4 Pi 

G p ,np(s, t, u, si, s 2 , s 3 , s 4 ) = T 2 (l + e) V -j- 

i=o t 




Figure 1: The (a) planar and (b) non-planar double box. 

We would like to note at this point that our method is not a priori a numerical 
method. In principle, the functions obtained after step 3 of the algorithm can be 
integrated analytically. However, as the parameter integrals for the coefficient of the 
l/e l pole are (N — I — Z)-dimensional, and the functions to integrate are in general 
highly non-trivial, only the integrals for large I, i.e. the highest pole coefficients, can 
be evaluated automatically by standard algebraic integration packages. For example, 
one obtains the following analytical result for the leading and subleading poles of the 
planar topology with s 3 7^ 0, s 4 7^ : 



G P (M,u,0,0,s 3 ,s 4 ) = r2 ( 1 + e )(- s )" 26 ^{ ^4 

-log(-t) + ilog(-s 3 ) + ilog(-s 4 )]} + 0{\ 



e 3 



e 2 ' 



Although the analytic integration could be pushed further by writing algebraic inte- 
gration routines which are specialised to deal with the functions occurring in these 
integrals, we do not pursue the analytic approach further here. 
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-*,-«) = (2/3, 2/3, 2/3) 




(-Si, -S 2 , -S 3 , 


Si) 


(1,0,0,1) 


(0,1,0,1) 


(0,0,1,1) 


Pa 




-0.8437 





-0.8438 


P3 




-2.052 





-2.053 


P2 




-2.189 


-3.571 


-10.53 


Pi 




4.394 


-6.585 


-48.60 


Po 




35.64 


11.83 


-140.7 




(-«, 


-t,-u) = (1/2, 1/3, 5/6) 




{ — s lj ~ s 2, —S3, 


-s 4 ) 


(2/3,0,0,1) 


(0,2/3,0,1) 


(0,0,2/3,1) 


P4 




-3.000 





-3.000 


P3 




-12.47 





-14.91 


P2 




-26.28 


-7.827 


-73.79 


Pi 




-19.04 


-18.90 


-325.3 


Po 




90.68 


16.18 


-1151. 



Table 1: Numerical results for the planar double box with two off-shell legs 

3.3 Double Box graphs for Bhabha scattering 

Bhabha scattering is a paradigm process of QED. Its calculation at the two- loop level, 
apart from the theoretical challenge alone, has a strong motivation by the fact that 
Bhabha scattering serves as a luminosity monitor for e + e~ colliders. 

The virtual two-loop QED corrections to Bhabha scattering have been calculated 
analytically by Bern, Dixon and Ghinculov [S] , in the limit of vanishing electron mass. 
Here we give numerical results for the three 2-loop box master integrals entering 
these corrections with massive internal lines. The numerical points are calculated for 
nonphysical kinematics, s, t < 0, but they can serve as a strong check for a future 
analytic calculation which includes the fermion masses. 

The three master topologies G/ =a ,fe,c are shown in Fig. El The corresponding func- 
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-*,-«) = (2/3, 2/3, 2/3) 




(-Si, -S 2 , -S 3 , 


Si) 


(1,1,0,0) 


(0,1,0,1) 


(0,0,1,1) 


Pa 




-2.812 


-0.5625 





P3 




-3.192 


-0.0923 


2.189 


P2 




29.82 


5.481 


1.634 


Pi 




139.2 


35.40 


-17.98 


Po 




365.3 


170.6 


-60.33 




(-«, 


-t,-u) = (1/2, 1/3, 5/6) 




{ — s lj ~ s 2, —S3, 


-s 4 ) 


(2/3,1,0,0) 


(0,2/3,0,1) 


(0,0,2/3,1) 


P4 




-6.700 


-0.600 





P3 




-15.97 


1.377 


4.502 


P2 




38.24 


11.79 


5.870 


Pi 




305.1 


57.33 


-36.18 


Po 




1080. 


267.0 


-165.5 



Table 2: Numerical results for the non-planar double box with two off-shell legs 

tions U and T are given by equation (jlj) for G a and G c and by © for G b : 

G a {s,t,u,m 2 ,M 2 ) = G P {s,t,u,m 2 ,M 2 ,M 2 ,m 2 ;0,m 2 ,M 2 ,0,m 2 ,M 2 ,0) 
G b {s,t,u,m 2 ,M 2 ) = G NP {s,t,u,m 2 ,M 2 ,M 2 ,m 2 ;0,m 2 ,M 2 ,0,m 2 ,M 2 ,0) 
G c {s,t,u,m 2 ,M 2 ) = Gp{s,t,u,m 2 ,m 2 ,M 2 ,M 2 ;0,m 2 ,m 2 ,m 2 ,0,0,M 2 ) 

Results for two different numerical points are given in Table H3 An overall factor 
T 2 (l + e) has been extracted: 

G l=aAc (s, t, u, m 2 , M 2 ) = T 2 (l + e) ^ ^ . 

i=o e 

To be slightly more general, the second point covers the case where the two massive 
propagator lines flowing through the graphs have different masses. 
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^wwvw\^ 



Kaaaa/vw^ 



Figure 2: The two-loop four point master topologies relevant for Bhabha scattering. 
The wavy lines are massless (photons) and the straight lines are massive scalars with 
external legs on-shell. We label the topologies from left to right by G a , Gb, G c . 



(s, -t,-u,m 2 ,M 2 ) 


(1/5,3/10,7/2,1,1) 


(5/3,4/3,5,1 


,3) 




G a 


Gb 


G c 


G a 


Gb 


G c 


P2 


-1.561 


-0.5255 


-1.152 


-0.08622 


-0.03483 


-0.05832 


Pi 


-5.335 


-0.2024 


-3.690 


-0.04195 


0.07556 


0.05389 


Po 


1.421 


3.606 


1.555 


0.7323 


0.1073 


0.6847 



Table 3: Results for the double box graphs for Bhabha scattering 
3.4 Propagators up to 5 loops 

Due to the high complexity of multi-loop QED and QCD calculations, only a few 
results at the three- and four loop level are known (for a review see e.g. [28J). These 
results rely to a large extent on the knowledge of multi-loop propagator functions. 

Massless propagator functions are very simple objects from the kinematical point 
of view, as they depend only on one single scale, s = p 2 , where p is the propagator 
momentum. Each graph is simply given by the scale to some power times a number 
which can be calculated analytically or numerically once and forever. 

We provide here some examples, shown in Figs. El to to demonstrate that our 
method can deal with different kinds of propagator topologies up to 5 loops, and 
that the treatment of UV sub divergences and higher order terms in the e-expansion 
is straightforward. 

For the planar and non-planar 3-loop ladder (see Fig. El), we get: 

G[3a] = (-s)- 2 ~ 3 T(2 + 3e) [20.74 + 93.71 e] (6) 
G[3b] = (-s)- 2 " 3 T(2 + 3e) [20.74 + 128.4 e] (7) 

For massless internal lines, UV sub divergences appear in the form of one-loop bubble 
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Figure 3: 3-loop ladder-type (a) planar and (b) non-planar propagator graph. 



insertions. The most convenient way to deal with them is to integrate them out 
analytically, which leads to a non-integer exponent of an internal propagator. An 
example is given in Fig.EJi : 



G[4a] = --(-s)- x - 4 T(l + 4e)Beta(l - e, 1 - e) 



20.74 + 86.57 e + 494.5 e' 



Note that in all three examples above the number 20 • £(5) = 20.73855. . . appears 
with the given precision 2 . Planar 4-loop ladder, Fig. 01) : 



(b) 



Figure 4: J^-loop propagators (a) with and (b) without UV subdivergence. 

G[4b] = -(-s)- 3 ~ 4 T(3 + 4e) [35.10 + 197.34 e] (9) 

The results for these propagator graphs agree with the analytical ones where avail- 
able (2ZI- Finally, we also calculated a 5- loop example, shown in Fig. 03 We obtain: 

G{5] = (-s)~ 4 - 5 T(4 + 5e) [40.53] (10) 




3.5 Massless on-shell planar 3-loop box 

In a very recent article [21] (see also [20j). Smirnov presented the analytical result 
for the massless on-shell planar 3-loop box TB(s,t), shown in Fig-El We have cross- 
checked the analytical result by calculating a numerical value for the point (s,t) = 
( — 1, —3), and obtain 

Going to higher precision is only a question of computer time, as the integrands are positive 
definite and bounded. 
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Figure 5: A 5-loop propagator graph 



Figure 6: Planar triple-box graph 



TB(-l,-3) = T(4 + 3e) 



r0.09874 0.7669 1.977 



0.7534 4.747 2.010 



+ 21.48 



which is in agreement within the 1% level with the analytical result. On a PC with a 
2 GHz Pentium IV processor, the full calculation took several weeks, but could have 
been speed up considerably by further parallelisation. However, disk space starts to 
become an issue here, as the size of the object files adds up to about 500 Mega Bytes. 



4 Discussion and Conclusions 

We have refined the formalism and automated program developed in pQ to isolate 
overlapping IR and UV divergences to be applicable to a very large class of Feyn- 
man graphs. The algorithm now is able to treat diagrams with arbitrary masses and 
propagator exponents. As discussed above, this also allows to compute Feynman 
integrals with non-trivial numerators. Hence we have formulated a constructive ap- 
proach, based on sector decomposition, to convert dimensionally regulated Feynman 
diagrams into a Laurent series in e, where the coefficients are given in terms of pa- 
rameter integrals. These parameter integrals can always be evaluated numerically in 
kinematic regions where the Mandelstam variables are negative, as one can show that 
the integrands are bounded and positive definite in that case. As examples, we calcu- 
lated numerically various 2-loop 4-point functions and a 3-loop 4-point function for 
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negative Mandelstam variables s, t, u. We also evaluated massless propagator graphs 
up to 5 loops. More precisely, we give numerical results for the following types of 
diagrams: 

- Massless two-loop 4-point functions with two off-shell external legs (planar as 
well as non-planar topologies). 

- The master two-loop box diagrams needed to calculate Bhabha scattering (with 
massive internal lines) at two loops. 

- Two-point functions at 3, 4 and 5 loops, among these a 4-loop graph containing 
an UV sub divergence, leading to non-integer propagator powers, as well as the 
order e terms for the 3-loop and 4-loop graphs. 

- The planar 3-loop massless 4-point function with on-shell legs. 

Integrals depending on a single scale, as for example massless 2-point functions, 
or massless 3-point functions with two on-shell legs, can be evaluated numerically if 
no analytical result is achievable, as they are just a number times the overall scale 
factor. We note that the algorithm can also be applied to vacuum graphs with UV 
sub divergences. In that case one has T = MJ2jLi x j m p an d any graph can be 
expressed by well-behaved integrals. In this special situation other very powerful 
methods are known |29| . 

For general multi-scale problems the situation is more delicate, as for physical kine- 
matics integrable singularities are present in a multi-dimensional integration space. 
The numerical integration in this case is a highly non-trivial task. For the one-loop 
case a combination of the sector decomposition algorithm and new numerical integra- 
tion methods were proposed in [TT] . where it has been exploited that sector decom- 
position also leads to more stable integral representations if no IR/UV singularities 
are present. A generalisation to the two-loop case is presently under investigation. 
Once a procedure is set up to deal with the sector integrals numerically for gen- 
eral kinematics, a numerical approach for multi-leg and multi-loop processes relevant 
at present and future high energy experiments will be feasible. In the meantime, 
nonetheless, the algorithm serves as an independent check for analytical results of 
multi-loop integrals. 
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